!
!     unsharp     convolution of the image by the gaussian mask
!     Copyright (C) 2004  Filip Hroch, Masaryk University, Brno, CZ
!
! $Id$
!

       Program Unsharp
!
! Experimental version 
!
! Shift frames.
! Using the output of the munimatch (with -t switch), eg. 'TRA'
! files with format:
! image1 image2 scale dx1, dx2 cos(fi) sin(fi) sin(fi) cos(fi)

     implicit none

! Parameters:

!    Version:
         Character(len=80),parameter :: version = 'Trashift v. 0.0.0'

! Initial parameters.
!    Output image name:
       Character(len=256) :: namekombi='unsharp.fits'

!    Minimum of values
       Real :: minvalue = 0.0

!    Maximum of values
       Real :: maxvalue = 65535.0

!-------

       Character(len=256) :: sb, sbtra
       Character(len=80) :: buf

       Logical :: simple,extend,anyf,exist
       Integer :: n, nalldate, maxalldate
       Integer :: i, j, istat, ostat ,bitpix, naxis, naxes(2) 
       Integer :: BlockSize,pcount,gcount

! image dimensions
       Integer ::  nx = -1, ny = -1 
! image
       Real, allocatable :: KOMBINE(:,:), IMAGE(:,:)
       real, allocatable :: mask(:,:)
       Real :: bscale, bzero, s, radius, ksum
       Real :: oneexp, exptime
       character(len=80) :: onedate
       character(len=80),allocatable :: alldate(:), alld(:)
       Real :: xx,yy,c,r(2),m(2,2),t1(2),t2(2)
       character(len=1000) :: sb1,sb2
       integer :: ix,jy, rdim

!------

       write(*,*) trim(version)

! init sum variables
       exptime = 0.0
       nalldate = 0
       maxalldate = 2
       allocate(alldate(maxalldate))

       ! nacteni jmena souboru a polomeru gaussovy masky 
       Read(*,*,ERR=666) sb, radius, c
       istat = 0

       ! initialisation of the mask
       rdim = max(2.0*radius, 3.0)
       allocate(mask(-rdim:rdim,-rdim:rdim))
       mask = 0
       ksum = 0.0
       do i = -rdim, rdim
          do j = -rdim, rdim
             mask(i,j) = exp(-(i**2 + j**2)/radius**2/2.0)
             ksum = ksum + mask(i,j)
          enddo
       enddo
!       mask = 0
       mask = -(c-1.0)*mask
       mask(0,0) = c*ksum -(c-1.0)*mask(0,0)
       

       ! input initialisation of the current file
       Call ftopen(25,sb,0,BlockSize,istat)
       Call ftghpr(25,2,simple,bitpix,naxis,naxes,pcount,gcount,extend,istat)
       if( istat /= 0 ) goto 666
       CALL ftgkye(25,'BSCALE',bscale,buf,istat)
       CALL ftgkye(25,'BZERO',bzero,buf,istat)
       !        BSCALE, BZERO not presented
       if( istat /= 0 ) then
          istat = 0
          bscale = 1.0
          bzero = 0.0
       endif
       call ftgkye(25,'EXPTIME',oneexp, buf, istat)
       if( istat /= 0 ) then
          oneexp = -1.0
          istat = 0
       endif
       call ftgkys(25,'DATE-OBS',onedate, buf, istat)
       if( istat /= 0 ) then
          onedate = ''
          istat = 0
       endif
       nx = naxes(1)
       ny = naxes(2)
       Allocate(   IMAGE(nx,ny) ) 
       IMAGE = 0.0
       Allocate( KOMBINE(nx,ny) )
       KOMBINE = 0.0
       call FTG2DE(25,1,minvalue,nx,nx,ny,IMAGE,anyf,istat)
       if( istat /= 0 ) goto 666
       call ftclos(25,istat)
       if( istat /= 0 ) goto 666

       ! compute convolution
       do i = rdim+1, nx - rdim - 1
          do j = rdim+1, ny - rdim - 1
             do ix = -rdim,rdim
                do jy = -rdim,rdim
                   KOMBINE(i,j) = KOMBINE(i,j) + mask(ix,jy)*IMAGE(ix+i,jy+j)
                enddo
             enddo
             if( abs(ksum) > epsilon(ksum) ) KOMBINE(i,j) = KOMBINE(i,j)/ksum
          enddo
       enddo
       deallocate(mask)

! form the output image
    ostat = 0
    Inquire(file=namekombi,exist=exist)
    if( exist ) then
       open(26,file=namekombi)
       close(26,status='DELETE')
    endif
    Call ftinit(26,namekombi,1,ostat)
       naxis = 2
       naxes = (/nx, ny/)
       bitpix = 32
       Call ftphps(26,bitpix,naxis,naxes,ostat)
!    call ftcopy(25,26,1,ostat)
    call ftphis(26,'Shifted image.',ostat)
    if( ostat /= 0 ) goto 666
    
    Call ftp2de(26,1,nx,nx,ny,KOMBINE,ostat)
    
    CALL ftclos(26,ostat)
    if( ostat /= 0 ) goto 666

! print output info
    write(*,*) 'Output image: ',trim(namekombi)
    write(*,*) 'Dimension: ',nx,'x',ny
    !       write(*,*) 'Bye.'
    stop 'Bye.'
    
666 continue 

    if( istat +  ostat /= 0 )then
       write(*,*) 'The following error(s) has been occured:'
       do
          Call ftgmsg(buf)
          if( buf == ' ' ) exit
          write(*,*) trim(buf)
       enddo
    endif
    
    stop 'Crached...'
    
  End       

